home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / ELECTRIC / DSPICE0S.ZIP / bjt.c < prev    next >
C/C++ Source or Header  |  1992-11-22  |  42KB  |  1,193 lines
  1. /* bjt.f -- translated by f2c (version of 3 February 1990  3:36:42).
  2.    You must link the resulting object file with the libraries:
  3.     -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #include "f2c.h"
  7.  
  8. /* Common Block Declarations */
  9.  
  10. struct {
  11.     integer ielmnt, isbckt, nsbckt, iunsat, nunsat, itemps, numtem, isens, 
  12.         nsens, ifour, nfour, ifield, icode, idelim, icolum, insize, 
  13.         junode, lsbkpt, numbkp, iorder, jmnode, iur, iuc, ilc, ilr, 
  14.         numoff, isr, nmoffc, iseq, iseq1, neqn, nodevs, ndiag, iswap, 
  15.         iequa, macins, lvnim1, lx0, lvn, lynl, lyu, lyl, lx1, lx2, lx3, 
  16.         lx4, lx5, lx6, lx7, ld0, ld1, ltd, imynl, imvn, lcvn, nsnod, 
  17.         nsmat, nsval, icnod, icmat, icval, loutpt, lpol, lzer, irswpf, 
  18.         irswpr, icswpf, icswpr, irpt, jcpt, irowno, jcolno, nttbr, nttar, 
  19.         lvntmp;
  20. } tabinf_;
  21.  
  22. #define tabinf_1 tabinf_
  23.  
  24. struct {
  25.     integer locate[50], jelcnt[50], nunods, ncnods, numnod, nstop, nut, nlt, 
  26.         nxtrm, ndist, ntlin, ibr, numvs, numalt, numcyc;
  27. } cirdat_;
  28.  
  29. #define cirdat_1 cirdat_
  30.  
  31. struct {
  32.     doublereal omega, time, delta, delold[7], ag[7], vt, xni, egfet, xmu, 
  33.         sfactr;
  34.     integer mode, modedc, icalc, initf, method, iord, maxord, noncon, iterno, 
  35.         itemno, nosolv, modac, ipiv, ivmflg, ipostp, iscrch, iofile;
  36. } status_;
  37.  
  38. #define status_1 status_
  39.  
  40. struct {
  41.     doublereal twopi, xlog2, xlog10, root2, rad, boltz, charge, ctok, gmin, 
  42.         reltol, abstol, vntol, trtol, chgtol, eps0, epssil, epsox, pivtol,
  43.          pivrel;
  44. } knstnt_;
  45.  
  46. #define knstnt_1 knstnt_
  47.  
  48. struct {
  49.     doublereal value[200000];
  50. } blank_;
  51.  
  52. #define blank_1 blank_
  53.  
  54. /*<       subroutine bjt >*/
  55. /* Subroutine */ int bjt_()
  56. {
  57.     /* System generated locals */
  58.     integer i_1;
  59.     doublereal d_1, d_2, d_3;
  60.  
  61.     /* Builtin functions */
  62.     double exp(), sqrt(), tan(), log();
  63.  
  64.     /* Local variables */
  65.     static doublereal cbcn, cben;
  66. #define cqbc ((doublereal *)&blank_1 + 11)
  67.     static doublereal area;
  68. #define cqbe ((doublereal *)&blank_1 + 9)
  69.     static doublereal gben, gbcn, fcpe, gccs, evbc, evbe, cdis;
  70.     static integer ioff;
  71.     static doublereal czbe, czbc;
  72. #define vbco ((doublereal *)&blank_1 + 1)
  73. #define vbeo ((doublereal *)&blank_1)
  74. #define cqcs ((doublereal *)&blank_1 + 13)
  75. #define gpio ((doublereal *)&blank_1 + 4)
  76. #define cqbx ((doublereal *)&blank_1 + 15)
  77. #define gmuo ((doublereal *)&blank_1 + 5)
  78.     static integer locv, locm, loct;
  79.     static doublereal csat, type, rbpr, rbpi, gcpr, gepr, oikr, xjrb, ctot, 
  80.         czbx, czcs, ovtf, xjtf, temp;
  81.     static integer ichk1;
  82.     static doublereal sarg, fcpc;
  83.     static integer locy, node1, node2, node3, node4, node5, node6, node7;
  84.     static doublereal capbc, capbe, ceqbc, ceqbe, geqcb;
  85. #define cexbc ((doublereal *)&blank_1 + 17)
  86.     static doublereal cchat, cbhat, capcs, ceqcs, evben, evbcn, ceqbx, denom, 
  87.         xfact, geqbx, argtf, capbx, sqarg, c2, c4, vcrit, q1, q2, f1, f2, 
  88.         f3, czbef2, czbcf2, czbxf2, cb, cc;
  89.     extern /* Subroutine */ int intgr8_();
  90.     static doublereal qb, gm, pe;
  91.     static integer icheck;
  92.     static doublereal go, td, tf, pc, delvbc;
  93. #define geqcbo ((doublereal *)&blank_1 + 18)
  94.     static doublereal delvbe, gx, dqbdve, dqbdvc, tr;
  95. #define nodplc ((integer *)&blank_1)
  96. #define cvalue ((complex *)&blank_1)
  97.     extern /* Subroutine */ int pnjlim_();
  98.     static doublereal cbc, ps, cbe, gbc, gbe;
  99. #define cbo ((doublereal *)&blank_1 + 3)
  100. #define cco ((doublereal *)&blank_1 + 2)
  101. #define qbc ((doublereal *)&blank_1 + 10)
  102. #define qbe ((doublereal *)&blank_1 + 8)
  103.     static doublereal bfm, vbc, arg;
  104.     static integer loc;
  105.     static doublereal vbe, brm;
  106. #define gmo ((doublereal *)&blank_1 + 6)
  107. #define goo ((doublereal *)&blank_1 + 7)
  108.     static doublereal ova;
  109. #define qcs ((doublereal *)&blank_1 + 12)
  110.     static doublereal ovb, oik;
  111. #define qbx ((doublereal *)&blank_1 + 14)
  112.     static doublereal vce, vtc;
  113. #define gxo ((doublereal *)&blank_1 + 16)
  114.     static doublereal vte, vbx, vcs, tol, gpi, gmu, vtn, cex, gex, xme, xmc, 
  115.         xms, xtf, geq, ceq, arg1, arg2, arg3;
  116.  
  117. /*<       implicit double precision (a-h,o-z) >*/
  118.  
  119. /*     this routine processes bjts for dc and transient analyses. */
  120.  
  121. /* spice version 2g.6  sccsid=tabinf 3/15/83 */
  122. /*<       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, >*/
  123. /*<      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, >*/
  124. /*<      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, >*/
  125. /*<      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, >*/
  126. /*<      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, >*/
  127. /*<      5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval, >*/
  128. /*<      6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt, >*/
  129. /*<      7   irowno,jcolno,nttbr,nttar,lvntmp >*/
  130. /* spice version 2g.6  sccsid=cirdat 3/15/83 */
  131. /*<       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop, >*/
  132. /*<      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc >*/
  133. /* spice version 2g.6  sccsid=status 3/15/83 */
  134. /*<       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, >*/
  135. /*<      1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon, >*/
  136. /*<      2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile >*/
  137. /* spice version 2g.6  sccsid=knstnt 3/15/83 */
  138. /*<       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok, >*/
  139. /*<      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox, >*/
  140. /*<      2   pivtol,pivrel >*/
  141. /* spice version 2g.6  sccsid=blank 3/15/83 */
  142. /*<       common /blank/ value(200000) >*/
  143. /*<       integer nodplc(64) >*/
  144. /*<       complex cvalue(32) >*/
  145. /*<       equivalence (value(1),nodplc(1),cvalue(1)) >*/
  146.  
  147.  
  148. /*<       dimension vbeo(1),vbco(1),cco(1),cbo(1),gpio(1),gmuo(1),gmo(1), >*/
  149. /*<      1   goo(1),qbe(1),cqbe(1),qbc(1),cqbc(1),qcs(1),cqcs(1),qbx(1), >*/
  150. /*<      2   cqbx(1),gxo(1),cexbc(1),geqcbo(1) >*/
  151. /*<       equivalence (vbeo(1),value(1)),(vbco(1),value(2)), >*/
  152. /*<      1   (cco(1),value(3)),(cbo(1),value(4)),(gpio(1),value(5)), >*/
  153. /*<      2   (gmuo(1),value(6)),(gmo(1),value(7)),(goo(1),value(8)), >*/
  154. /*<      3   (qbe(1),value(9)),(cqbe(1),value(10)),(qbc(1),value(11)), >*/
  155. /*<      4   (cqbc(1),value(12)),(qcs(1),value(13)),(cqcs(1),value(14)), >*/
  156. /*<      5   (qbx(1),value(15)),(cqbx(1),value(16)),(gxo(1),value(17)), >*/
  157. /*<      6   (cexbc(1),value(18)),(geqcbo(1),value(19)) >*/
  158.  
  159.  
  160. /*<       loc=locate(12) >*/
  161.     loc = cirdat_1.locate[11];
  162. /*<    10 if ((loc.eq.0).or.(nodplc(loc+36).ne.0)) return >*/
  163. L10:
  164.     if (loc == 0 || nodplc[loc + 35] != 0) {
  165.     return 0;
  166.     }
  167. /*<       locv=nodplc(loc+1) >*/
  168.     locv = nodplc[loc];
  169. /*<       node1=nodplc(loc+2) >*/
  170.     node1 = nodplc[loc + 1];
  171. /*<       node2=nodplc(loc+3) >*/
  172.     node2 = nodplc[loc + 2];
  173. /*<       node3=nodplc(loc+4) >*/
  174.     node3 = nodplc[loc + 3];
  175. /*<       node4=nodplc(loc+5) >*/
  176.     node4 = nodplc[loc + 4];
  177. /*<       node5=nodplc(loc+6) >*/
  178.     node5 = nodplc[loc + 5];
  179. /*<       node6=nodplc(loc+7) >*/
  180.     node6 = nodplc[loc + 6];
  181. /*<       node7=nodplc(loc+30) >*/
  182.     node7 = nodplc[loc + 29];
  183. /*<       locm=nodplc(loc+8) >*/
  184.     locm = nodplc[loc + 7];
  185. /*<       ioff=nodplc(loc+9) >*/
  186.     ioff = nodplc[loc + 8];
  187. /*<       type=nodplc(locm+2) >*/
  188.     type = (doublereal) nodplc[locm + 1];
  189. /*<       locm=nodplc(locm+1) >*/
  190.     locm = nodplc[locm];
  191. /*<       loct=nodplc(loc+22) >*/
  192.     loct = nodplc[loc + 21];
  193. /*<       gccs=0.0d0 >*/
  194.     gccs = 0.;
  195. /*<       ceqcs=0.0d0 >*/
  196.     ceqcs = 0.;
  197. /*<       geqbx=0.0d0 >*/
  198.     geqbx = 0.;
  199. /*<       ceqbx=0.0d0 >*/
  200.     ceqbx = 0.;
  201. /*<       geqcb=0.0d0 >*/
  202.     geqcb = 0.;
  203.  
  204. /*  dc model paramters */
  205.  
  206. /*<       area=value(locv+1) >*/
  207.     area = blank_1.value[locv];
  208. /*<       bfm=value(locm+2) >*/
  209.     bfm = blank_1.value[locm + 1];
  210. /*<       brm=value(locm+8) >*/
  211.     brm = blank_1.value[locm + 7];
  212. /*<       csat=value(locm+1)*area >*/
  213.     csat = blank_1.value[locm] * area;
  214. /*<       rbpr=value(locm+18)/area >*/
  215.     rbpr = blank_1.value[locm + 17] / area;
  216. /*<       rbpi=value(locm+16)/area-rbpr >*/
  217.     rbpi = blank_1.value[locm + 15] / area - rbpr;
  218. /*<       gcpr=value(locm+20)*area >*/
  219.     gcpr = blank_1.value[locm + 19] * area;
  220. /*<       gepr=value(locm+19)*area >*/
  221.     gepr = blank_1.value[locm + 18] * area;
  222. /*<       ova=value(locm+4) >*/
  223.     ova = blank_1.value[locm + 3];
  224. /*<       ovb=value(locm+10) >*/
  225.     ovb = blank_1.value[locm + 9];
  226. /*<       oik=value(locm+5)/area >*/
  227.     oik = blank_1.value[locm + 4] / area;
  228. /*<       c2=value(locm+6)*area >*/
  229.     c2 = blank_1.value[locm + 5] * area;
  230. /*<       vte=value(locm+7)*vt >*/
  231.     vte = blank_1.value[locm + 6] * status_1.vt;
  232. /*<       oikr=value(locm+11)/area >*/
  233.     oikr = blank_1.value[locm + 10] / area;
  234. /*<       c4=value(locm+12)*area >*/
  235.     c4 = blank_1.value[locm + 11] * area;
  236. /*<       vtc=value(locm+13)*vt >*/
  237.     vtc = blank_1.value[locm + 12] * status_1.vt;
  238. /*<       vcrit=value(locm+54) >*/
  239.     vcrit = blank_1.value[locm + 53];
  240. /*<       td=value(locm+28) >*/
  241.     td = blank_1.value[locm + 27];
  242. /*<       xjrb=value(locm+17)*area >*/
  243.     xjrb = blank_1.value[locm + 16] * area;
  244.  
  245. /*  initialization */
  246.  
  247. /*<       icheck=1 >*/
  248.     icheck = 1;
  249. /*<       go to (100,20,30,50,60,70),initf >*/
  250.     switch (status_1.initf) {
  251.     case 1:  goto L100;
  252.     case 2:  goto L20;
  253.     case 3:  goto L30;
  254.     case 4:  goto L50;
  255.     case 5:  goto L60;
  256.     case 6:  goto L70;
  257.     }
  258. /*<    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25 >*/
  259. L20:
  260.     if (status_1.mode != 1 || status_1.modedc != 2 || status_1.nosolv == 0) {
  261.     goto L25;
  262.     }
  263. /*<       vbe=type*value(locv+2) >*/
  264.     vbe = type * blank_1.value[locv + 1];
  265. /*<       vce=type*value(locv+3) >*/
  266.     vce = type * blank_1.value[locv + 2];
  267. /*<       vbc=vbe-vce >*/
  268.     vbc = vbe - vce;
  269. /*<       vbx=vbc >*/
  270.     vbx = vbc;
  271. /*<       vcs=0.0d0 >*/
  272.     vcs = 0.;
  273. /*<       go to 300 >*/
  274.     goto L300;
  275. /*<    25 if(ioff.ne.0) go to 40 >*/
  276. L25:
  277.     if (ioff != 0) {
  278.     goto L40;
  279.     }
  280. /*<       vbe=vcrit >*/
  281.     vbe = vcrit;
  282. /*<       vbc=0.0d0 >*/
  283.     vbc = 0.;
  284. /*<       go to 300 >*/
  285.     goto L300;
  286. /*<    30 if (ioff.eq.0) go to 100 >*/
  287. L30:
  288.     if (ioff == 0) {
  289.     goto L100;
  290.     }
  291. /*<    40 vbe=0.0d0 >*/
  292. L40:
  293.     vbe = 0.;
  294. /*<       vbc=0.0d0 >*/
  295.     vbc = 0.;
  296. /*<       go to 300 >*/
  297.     goto L300;
  298. /*<    50 vbe=vbeo(lx0+loct) >*/
  299. L50:
  300.     vbe = vbeo[tabinf_1.lx0 + loct - 1];
  301. /*<       vbc=vbco(lx0+loct) >*/
  302.     vbc = vbco[tabinf_1.lx0 + loct - 1];
  303. /*<       vbx=type*(value(lvnim1+node2)-value(lvnim1+node4)) >*/
  304.     vbx = type * (blank_1.value[tabinf_1.lvnim1 + node2 - 1] - blank_1.value[
  305.         tabinf_1.lvnim1 + node4 - 1]);
  306. /*<       vcs=type*(value(lvnim1+node7)-value(lvnim1+node4)) >*/
  307.     vcs = type * (blank_1.value[tabinf_1.lvnim1 + node7 - 1] - blank_1.value[
  308.         tabinf_1.lvnim1 + node4 - 1]);
  309. /*<       go to 300 >*/
  310.     goto L300;
  311. /*<    60 vbe=vbeo(lx1+loct) >*/
  312. L60:
  313.     vbe = vbeo[tabinf_1.lx1 + loct - 1];
  314. /*<       vbc=vbco(lx1+loct) >*/
  315.     vbc = vbco[tabinf_1.lx1 + loct - 1];
  316. /*<       vbx=type*(value(lvnim1+node2)-value(lvnim1+node4)) >*/
  317.     vbx = type * (blank_1.value[tabinf_1.lvnim1 + node2 - 1] - blank_1.value[
  318.         tabinf_1.lvnim1 + node4 - 1]);
  319. /*<       vcs=type*(value(lvnim1+node7)-value(lvnim1+node4)) >*/
  320.     vcs = type * (blank_1.value[tabinf_1.lvnim1 + node7 - 1] - blank_1.value[
  321.         tabinf_1.lvnim1 + node4 - 1]);
  322. /*<       if(mode.ne.2.or.nosolv.eq.0) go to 300 >*/
  323.     if (status_1.mode != 2 || status_1.nosolv == 0) {
  324.     goto L300;
  325.     }
  326. /*<       vbx=type*(value(locv+2)-value(locv+3)) >*/
  327.     vbx = type * (blank_1.value[locv + 1] - blank_1.value[locv + 2]);
  328. /*<       vcs=0.0d0 >*/
  329.     vcs = 0.;
  330. /*<       go to 300 >*/
  331.     goto L300;
  332. /*<    70 xfact=delta/delold(2) >*/
  333. L70:
  334.     xfact = status_1.delta / status_1.delold[1];
  335. /*<       vbeo(lx0+loct)=vbeo(lx1+loct) >*/
  336.     vbeo[tabinf_1.lx0 + loct - 1] = vbeo[tabinf_1.lx1 + loct - 1];
  337. /*<       vbe=(1.0d0+xfact)*vbeo(lx1+loct)-xfact*vbeo(lx2+loct) >*/
  338.     vbe = (xfact + 1.) * vbeo[tabinf_1.lx1 + loct - 1] - xfact * vbeo[
  339.         tabinf_1.lx2 + loct - 1];
  340. /*<       vbco(lx0+loct)=vbco(lx1+loct) >*/
  341.     vbco[tabinf_1.lx0 + loct - 1] = vbco[tabinf_1.lx1 + loct - 1];
  342. /*<       vbc=(1.0d0+xfact)*vbco(lx1+loct)-xfact*vbco(lx2+loct) >*/
  343.     vbc = (xfact + 1.) * vbco[tabinf_1.lx1 + loct - 1] - xfact * vbco[
  344.         tabinf_1.lx2 + loct - 1];
  345. /*<       cco(lx0+loct)=cco(lx1+loct) >*/
  346.     cco[tabinf_1.lx0 + loct - 1] = cco[tabinf_1.lx1 + loct - 1];
  347. /*<       cbo(lx0+loct)=cbo(lx1+loct) >*/
  348.     cbo[tabinf_1.lx0 + loct - 1] = cbo[tabinf_1.lx1 + loct - 1];
  349. /*<       gpio(lx0+loct)=gpio(lx1+loct) >*/
  350.     gpio[tabinf_1.lx0 + loct - 1] = gpio[tabinf_1.lx1 + loct - 1];
  351. /*<       gmuo(lx0+loct)=gmuo(lx1+loct) >*/
  352.     gmuo[tabinf_1.lx0 + loct - 1] = gmuo[tabinf_1.lx1 + loct - 1];
  353. /*<       gmo(lx0+loct)=gmo(lx1+loct) >*/
  354.     gmo[tabinf_1.lx0 + loct - 1] = gmo[tabinf_1.lx1 + loct - 1];
  355. /*<       goo(lx0+loct)=goo(lx1+loct) >*/
  356.     goo[tabinf_1.lx0 + loct - 1] = goo[tabinf_1.lx1 + loct - 1];
  357. /*<       gxo(lx0+loct)=gxo(lx1+loct) >*/
  358.     gxo[tabinf_1.lx0 + loct - 1] = gxo[tabinf_1.lx1 + loct - 1];
  359. /*<       go to 110 >*/
  360.     goto L110;
  361.  
  362. /*  compute new nonlinear branch voltages */
  363.  
  364. /*<   100 vbe=type*(value(lvnim1+node5)-value(lvnim1+node6)) >*/
  365. L100:
  366.     vbe = type * (blank_1.value[tabinf_1.lvnim1 + node5 - 1] - blank_1.value[
  367.         tabinf_1.lvnim1 + node6 - 1]);
  368. /*<       vbc=type*(value(lvnim1+node5)-value(lvnim1+node4)) >*/
  369.     vbc = type * (blank_1.value[tabinf_1.lvnim1 + node5 - 1] - blank_1.value[
  370.         tabinf_1.lvnim1 + node4 - 1]);
  371. /*<   110 delvbe=vbe-vbeo(lx0+loct) >*/
  372. L110:
  373.     delvbe = vbe - vbeo[tabinf_1.lx0 + loct - 1];
  374. /*<       delvbc=vbc-vbco(lx0+loct) >*/
  375.     delvbc = vbc - vbco[tabinf_1.lx0 + loct - 1];
  376. /*<       vbx=type*(value(lvnim1+node2)-value(lvnim1+node4)) >*/
  377.     vbx = type * (blank_1.value[tabinf_1.lvnim1 + node2 - 1] - blank_1.value[
  378.         tabinf_1.lvnim1 + node4 - 1]);
  379. /*<       vcs=type*(value(lvnim1+node7)-value(lvnim1+node4)) >*/
  380.     vcs = type * (blank_1.value[tabinf_1.lvnim1 + node7 - 1] - blank_1.value[
  381.         tabinf_1.lvnim1 + node4 - 1]);
  382. /*<       cchat=cco(lx0+loct)+(gmo(lx0+loct)+goo(lx0+loct))*delvbe >*/
  383. /*<      1   -(goo(lx0+loct)+gmuo(lx0+loct))*delvbc >*/
  384.     cchat = cco[tabinf_1.lx0 + loct - 1] + (gmo[tabinf_1.lx0 + loct - 1] + 
  385.         goo[tabinf_1.lx0 + loct - 1]) * delvbe - (goo[tabinf_1.lx0 + loct 
  386.         - 1] + gmuo[tabinf_1.lx0 + loct - 1]) * delvbc;
  387. /*<       cbhat=cbo(lx0+loct)+gpio(lx0+loct)*delvbe+gmuo(lx0+loct)*delvbc >*/
  388.     cbhat = cbo[tabinf_1.lx0 + loct - 1] + gpio[tabinf_1.lx0 + loct - 1] * 
  389.         delvbe + gmuo[tabinf_1.lx0 + loct - 1] * delvbc;
  390.  
  391. /*   bypass if solution has not changed */
  392.  
  393. /*<       if (initf.eq.6) go to 200 >*/
  394.     if (status_1.initf == 6) {
  395.     goto L200;
  396.     }
  397. /*<       tol=reltol*dmax1(dabs(vbe),dabs(vbeo(lx0+loct)))+vntol >*/
  398. /* Computing MAX */
  399.     d_2 = abs(vbe), d_3 = (d_1 = vbeo[tabinf_1.lx0 + loct - 1], abs(d_1));
  400.     tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.vntol;
  401. /*<       if (dabs(delvbe).ge.tol) go to 200 >*/
  402.     if (abs(delvbe) >= tol) {
  403.     goto L200;
  404.     }
  405. /*<       tol=reltol*dmax1(dabs(vbc),dabs(vbco(lx0+loct)))+vntol >*/
  406. /* Computing MAX */
  407.     d_2 = abs(vbc), d_3 = (d_1 = vbco[tabinf_1.lx0 + loct - 1], abs(d_1));
  408.     tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.vntol;
  409. /*<       if (dabs(delvbc).ge.tol) go to 200 >*/
  410.     if (abs(delvbc) >= tol) {
  411.     goto L200;
  412.     }
  413. /*<       tol=reltol*dmax1(dabs(cchat),dabs(cco(lx0+loct)))+abstol >*/
  414. /* Computing MAX */
  415.     d_2 = abs(cchat), d_3 = (d_1 = cco[tabinf_1.lx0 + loct - 1], abs(d_1));
  416.     tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.abstol;
  417. /*<       if (dabs(cchat-cco(lx0+loct)).ge.tol) go to 200 >*/
  418.     if ((d_1 = cchat - cco[tabinf_1.lx0 + loct - 1], abs(d_1)) >= tol) {
  419.     goto L200;
  420.     }
  421. /*<       tol=reltol*dmax1(dabs(cbhat),dabs(cbo(lx0+loct)))+abstol >*/
  422. /* Computing MAX */
  423.     d_2 = abs(cbhat), d_3 = (d_1 = cbo[tabinf_1.lx0 + loct - 1], abs(d_1));
  424.     tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.abstol;
  425. /*<       if (dabs(cbhat-cbo(lx0+loct)).ge.tol) go to 200 >*/
  426.     if ((d_1 = cbhat - cbo[tabinf_1.lx0 + loct - 1], abs(d_1)) >= tol) {
  427.     goto L200;
  428.     }
  429. /*<       vbe=vbeo(lx0+loct) >*/
  430.     vbe = vbeo[tabinf_1.lx0 + loct - 1];
  431. /*<       vbc=vbco(lx0+loct) >*/
  432.     vbc = vbco[tabinf_1.lx0 + loct - 1];
  433. /*<       cc=cco(lx0+loct) >*/
  434.     cc = cco[tabinf_1.lx0 + loct - 1];
  435. /*<       cb=cbo(lx0+loct) >*/
  436.     cb = cbo[tabinf_1.lx0 + loct - 1];
  437. /*<       gpi=gpio(lx0+loct) >*/
  438.     gpi = gpio[tabinf_1.lx0 + loct - 1];
  439. /*<       gmu=gmuo(lx0+loct) >*/
  440.     gmu = gmuo[tabinf_1.lx0 + loct - 1];
  441. /*<       gm=gmo(lx0+loct) >*/
  442.     gm = gmo[tabinf_1.lx0 + loct - 1];
  443. /*<       go=goo(lx0+loct) >*/
  444.     go = goo[tabinf_1.lx0 + loct - 1];
  445. /*<       gx=gxo(lx0+loct) >*/
  446.     gx = gxo[tabinf_1.lx0 + loct - 1];
  447. /*<       geqcb=geqcbo(lx0+loct) >*/
  448.     geqcb = geqcbo[tabinf_1.lx0 + loct - 1];
  449. /*<       if (mode.ne.1) go to 800 >*/
  450.     if (status_1.mode != 1) {
  451.     goto L800;
  452.     }
  453. /*<       go to 900 >*/
  454.     goto L900;
  455.  
  456. /*  limit nonlinear branch voltages */
  457.  
  458. /*<   200 ichk1=1 >*/
  459. L200:
  460.     ichk1 = 1;
  461. /*<       call pnjlim(vbe,vbeo(lx0+loct),vt,vcrit,icheck) >*/
  462.     pnjlim_(&vbe, &vbeo[tabinf_1.lx0 + loct - 1], &status_1.vt, &vcrit, &
  463.         icheck);
  464. /*<       call pnjlim(vbc,vbco(lx0+loct),vt,vcrit,ichk1) >*/
  465.     pnjlim_(&vbc, &vbco[tabinf_1.lx0 + loct - 1], &status_1.vt, &vcrit, &
  466.         ichk1);
  467. /*<       if (ichk1.eq.1) icheck=1 >*/
  468.     if (ichk1 == 1) {
  469.     icheck = 1;
  470.     }
  471.  
  472. /*  determine dc current and derivitives */
  473.  
  474. /*<   300 vtn=vt*value(locm+3) >*/
  475. L300:
  476.     vtn = status_1.vt * blank_1.value[locm + 2];
  477. /*<       if(vbe.le.-5.0d0*vtn) go to 320 >*/
  478.     if (vbe <= vtn * -5.) {
  479.     goto L320;
  480.     }
  481. /*<       evbe=dexp(vbe/vtn) >*/
  482.     evbe = exp(vbe / vtn);
  483. /*<       cbe=csat*(evbe-1.0d0)+gmin*vbe >*/
  484.     cbe = csat * (evbe - 1.) + knstnt_1.gmin * vbe;
  485. /*<       gbe=csat*evbe/vtn+gmin >*/
  486.     gbe = csat * evbe / vtn + knstnt_1.gmin;
  487. /*<       if (c2.ne.0.0d0) go to 310 >*/
  488.     if (c2 != 0.) {
  489.     goto L310;
  490.     }
  491. /*<       cben=0.0d0 >*/
  492.     cben = 0.;
  493. /*<       gben=0.0d0 >*/
  494.     gben = 0.;
  495. /*<       go to 350 >*/
  496.     goto L350;
  497. /*<   310 evben=dexp(vbe/vte) >*/
  498. L310:
  499.     evben = exp(vbe / vte);
  500. /*<       cben=c2*(evben-1.0d0) >*/
  501.     cben = c2 * (evben - 1.);
  502. /*<       gben=c2*evben/vte >*/
  503.     gben = c2 * evben / vte;
  504. /*<       go to 350 >*/
  505.     goto L350;
  506. /*<   320 gbe=-csat/vbe+gmin >*/
  507. L320:
  508.     gbe = -csat / vbe + knstnt_1.gmin;
  509. /*<       cbe=gbe*vbe >*/
  510.     cbe = gbe * vbe;
  511. /*<       gben=-c2/vbe >*/
  512.     gben = -c2 / vbe;
  513. /*<       cben=gben*vbe >*/
  514.     cben = gben * vbe;
  515. /*<   350 vtn=vt*value(locm+9) >*/
  516. L350:
  517.     vtn = status_1.vt * blank_1.value[locm + 8];
  518. /*<       if(vbc.le.-5.0d0*vtn) go to 370 >*/
  519.     if (vbc <= vtn * -5.) {
  520.     goto L370;
  521.     }
  522. /*<       evbc=dexp(vbc/vtn) >*/
  523.     evbc = exp(vbc / vtn);
  524. /*<       cbc=csat*(evbc-1.0d0)+gmin*vbc >*/
  525.     cbc = csat * (evbc - 1.) + knstnt_1.gmin * vbc;
  526. /*<       gbc=csat*evbc/vtn+gmin >*/
  527.     gbc = csat * evbc / vtn + knstnt_1.gmin;
  528. /*<       if (c4.ne.0.0d0) go to 360 >*/
  529.     if (c4 != 0.) {
  530.     goto L360;
  531.     }
  532. /*<       cbcn=0.0d0 >*/
  533.     cbcn = 0.;
  534. /*<       gbcn=0.0d0 >*/
  535.     gbcn = 0.;
  536. /*<       go to 400 >*/
  537.     goto L400;
  538. /*<   360 evbcn=dexp(vbc/vtc) >*/
  539. L360:
  540.     evbcn = exp(vbc / vtc);
  541. /*<       cbcn=c4*(evbcn-1.0d0) >*/
  542.     cbcn = c4 * (evbcn - 1.);
  543. /*<       gbcn=c4*evbcn/vtc >*/
  544.     gbcn = c4 * evbcn / vtc;
  545. /*<       go to 400 >*/
  546.     goto L400;
  547. /*<   370 gbc=-csat/vbc+gmin >*/
  548. L370:
  549.     gbc = -csat / vbc + knstnt_1.gmin;
  550. /*<       cbc=gbc*vbc >*/
  551.     cbc = gbc * vbc;
  552. /*<       gbcn=-c4/vbc >*/
  553.     gbcn = -c4 / vbc;
  554. /*<       cbcn=gbcn*vbc >*/
  555.     cbcn = gbcn * vbc;
  556.  
  557. /*  determine base charge terms */
  558.  
  559. /*<   400 q1=1.0d0/(1.0d0-ova*vbc-ovb*vbe) >*/
  560. L400:
  561.     q1 = 1. / (1. - ova * vbc - ovb * vbe);
  562. /*<       if (oik.ne.0.0d0) go to 405 >*/
  563.     if (oik != 0.) {
  564.     goto L405;
  565.     }
  566. /*<       if (oikr.ne.0.0d0) go to 405 >*/
  567.     if (oikr != 0.) {
  568.     goto L405;
  569.     }
  570. /*<       qb=q1 >*/
  571.     qb = q1;
  572. /*<       dqbdve=q1*qb*ovb >*/
  573.     dqbdve = q1 * qb * ovb;
  574. /*<       dqbdvc=q1*qb*ova >*/
  575.     dqbdvc = q1 * qb * ova;
  576. /*<       go to 410 >*/
  577.     goto L410;
  578. /*<   405 q2=oik*cbe+oikr*cbc >*/
  579. L405:
  580.     q2 = oik * cbe + oikr * cbc;
  581. /*<       arg=dmax1(0.0d0,1.0d0+4.0d0*q2) >*/
  582. /* Computing MAX */
  583.     d_1 = 0., d_2 = q2 * 4. + 1.;
  584.     arg = max(d_2,d_1);
  585. /*<       sqarg=1.0d0 >*/
  586.     sqarg = 1.;
  587. /*<       if(arg.ne.0.0d0) sqarg=dsqrt(arg) >*/
  588.     if (arg != 0.) {
  589.     sqarg = sqrt(arg);
  590.     }
  591. /*<       qb=q1*(1.0d0+sqarg)/2.0d0 >*/
  592.     qb = q1 * (sqarg + 1.) / 2.;
  593. /*<       dqbdve=q1*(qb*ovb+oik*gbe/sqarg) >*/
  594.     dqbdve = q1 * (qb * ovb + oik * gbe / sqarg);
  595. /*<       dqbdvc=q1*(qb*ova+oikr*gbc/sqarg) >*/
  596.     dqbdvc = q1 * (qb * ova + oikr * gbc / sqarg);
  597.  
  598. /*  weil's approx. for excess phase applied with backward- */
  599. /*  euler integration */
  600.  
  601. /*<   410 cc=0.0d0 >*/
  602. L410:
  603.     cc = 0.;
  604. /*<       cex=cbe >*/
  605.     cex = cbe;
  606. /*<       gex=gbe >*/
  607.     gex = gbe;
  608. /*<       if(mode.eq.1) go to 420 >*/
  609.     if (status_1.mode == 1) {
  610.     goto L420;
  611.     }
  612. /*<       if(td.eq.0.0d0) go to 420 >*/
  613.     if (td == 0.) {
  614.     goto L420;
  615.     }
  616. /*<       arg1=delta/td >*/
  617.     arg1 = status_1.delta / td;
  618. /*<       arg2=3.0d0*arg1 >*/
  619.     arg2 = arg1 * 3.;
  620. /*<       arg1=arg2*arg1 >*/
  621.     arg1 = arg2 * arg1;
  622. /*<       denom=1.0d0+arg1+arg2 >*/
  623.     denom = arg1 + 1. + arg2;
  624. /*<       arg3=arg1/denom >*/
  625.     arg3 = arg1 / denom;
  626. /*<       if(initf.ne.5) go to 411 >*/
  627.     if (status_1.initf != 5) {
  628.     goto L411;
  629.     }
  630. /*<       cexbc(lx1+loct)=cbe/qb >*/
  631.     cexbc[tabinf_1.lx1 + loct - 1] = cbe / qb;
  632. /*<       cexbc(lx2+loct)=cexbc(lx1+loct) >*/
  633.     cexbc[tabinf_1.lx2 + loct - 1] = cexbc[tabinf_1.lx1 + loct - 1];
  634. /*<   411 cc=(cexbc(lx1+loct)*(1.0d0+delta/delold(2)+arg2) >*/
  635. /*<      1  -cexbc(lx2+loct)*delta/delold(2))/denom >*/
  636. L411:
  637.     cc = (cexbc[tabinf_1.lx1 + loct - 1] * (status_1.delta / status_1.delold[
  638.         1] + 1. + arg2) - cexbc[tabinf_1.lx2 + loct - 1] * status_1.delta 
  639.         / status_1.delold[1]) / denom;
  640. /*<       cex=cbe*arg3 >*/
  641.     cex = cbe * arg3;
  642. /*<       gex=gbe*arg3 >*/
  643.     gex = gbe * arg3;
  644. /*<       cexbc(lx0+loct)=cc+cex/qb >*/
  645.     cexbc[tabinf_1.lx0 + loct - 1] = cc + cex / qb;
  646.  
  647. /*  determine dc incremental conductances */
  648.  
  649. /*<   420 cc=cc+(cex-cbc)/qb-cbc/brm-cbcn >*/
  650. L420:
  651.     cc = cc + (cex - cbc) / qb - cbc / brm - cbcn;
  652. /*<       cb=cbe/bfm+cben+cbc/brm+cbcn >*/
  653.     cb = cbe / bfm + cben + cbc / brm + cbcn;
  654. /*<       gx=rbpr+rbpi/qb >*/
  655.     gx = rbpr + rbpi / qb;
  656. /*<       if(xjrb.eq.0.0d0) go to 430 >*/
  657.     if (xjrb == 0.) {
  658.     goto L430;
  659.     }
  660. /*<       arg1=dmax1(cb/xjrb,1.0d-9) >*/
  661. /* Computing MAX */
  662.     d_1 = cb / xjrb;
  663.     arg1 = max(1e-9,d_1);
  664. /*<       arg2=(-1.0d0+dsqrt(1.0d0+14.59025d0*arg1))/2.4317d0/dsqrt(arg1) >*/
  665.     arg2 = (sqrt(arg1 * 14.59025 + 1.) - 1.) / 2.4317 / sqrt(arg1);
  666. /*<       arg1=dtan(arg2) >*/
  667.     arg1 = tan(arg2);
  668. /*<       gx=rbpr+3.0d0*rbpi*(arg1-arg2)/arg2/arg1/arg1 >*/
  669.     gx = rbpr + rbpi * 3. * (arg1 - arg2) / arg2 / arg1 / arg1;
  670. /*<   430 if(gx.ne.0.0d0) gx=1.0d0/gx >*/
  671. L430:
  672.     if (gx != 0.) {
  673.     gx = 1. / gx;
  674.     }
  675. /*<       gpi=gbe/bfm+gben >*/
  676.     gpi = gbe / bfm + gben;
  677. /*<       gmu=gbc/brm+gbcn >*/
  678.     gmu = gbc / brm + gbcn;
  679. /*<       go=(gbc+(cex-cbc)*dqbdvc/qb)/qb >*/
  680.     go = (gbc + (cex - cbc) * dqbdvc / qb) / qb;
  681. /*<       gm=(gex-(cex-cbc)*dqbdve/qb)/qb-go >*/
  682.     gm = (gex - (cex - cbc) * dqbdve / qb) / qb - go;
  683. /*<       if (mode.ne.1) go to 500 >*/
  684.     if (status_1.mode != 1) {
  685.     goto L500;
  686.     }
  687. /*<       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500 >*/
  688.     if (status_1.modedc == 2 && status_1.nosolv != 0) {
  689.     goto L500;
  690.     }
  691. /*<       if (initf.eq.4) go to 500 >*/
  692.     if (status_1.initf == 4) {
  693.     goto L500;
  694.     }
  695. /*<       go to 700 >*/
  696.     goto L700;
  697.  
  698. /*  charge storage elements */
  699.  
  700. /*<   500 tf=value(locm+24) >*/
  701. L500:
  702.     tf = blank_1.value[locm + 23];
  703. /*<       tr=value(locm+33) >*/
  704.     tr = blank_1.value[locm + 32];
  705. /*<       czbe=value(locm+21)*area >*/
  706.     czbe = blank_1.value[locm + 20] * area;
  707. /*<       pe=value(locm+22) >*/
  708.     pe = blank_1.value[locm + 21];
  709. /*<       xme=value(locm+23) >*/
  710.     xme = blank_1.value[locm + 22];
  711. /*<       cdis=value(locm+32) >*/
  712.     cdis = blank_1.value[locm + 31];
  713. /*<       ctot=value(locm+29)*area >*/
  714.     ctot = blank_1.value[locm + 28] * area;
  715. /*<       czbc=ctot*cdis >*/
  716.     czbc = ctot * cdis;
  717. /*<       czbx=ctot-czbc >*/
  718.     czbx = ctot - czbc;
  719. /*<       pc=value(locm+30) >*/
  720.     pc = blank_1.value[locm + 29];
  721. /*<       xmc=value(locm+31) >*/
  722.     xmc = blank_1.value[locm + 30];
  723. /*<       fcpe=value(locm+46) >*/
  724.     fcpe = blank_1.value[locm + 45];
  725. /*<       czcs=value(locm+38)*area >*/
  726.     czcs = blank_1.value[locm + 37] * area;
  727. /*<       ps=value(locm+39) >*/
  728.     ps = blank_1.value[locm + 38];
  729. /*<       xms=value(locm+40) >*/
  730.     xms = blank_1.value[locm + 39];
  731. /*<       xtf=value(locm+25) >*/
  732.     xtf = blank_1.value[locm + 24];
  733. /*<       ovtf=value(locm+26) >*/
  734.     ovtf = blank_1.value[locm + 25];
  735. /*<       xjtf=value(locm+27)*area >*/
  736.     xjtf = blank_1.value[locm + 26] * area;
  737. /*<       if(tf.eq.0.0d0) go to 505 >*/
  738.     if (tf == 0.) {
  739.     goto L505;
  740.     }
  741. /*<       if(vbe.le.0.0d0) go to 505 >*/
  742.     if (vbe <= 0.) {
  743.     goto L505;
  744.     }
  745. /*<       argtf=0.0d0 >*/
  746.     argtf = 0.;
  747. /*<       arg2=0.0d0 >*/
  748.     arg2 = 0.;
  749. /*<       arg3=0.0d0 >*/
  750.     arg3 = 0.;
  751. /*<       if(xtf.eq.0.0d0) go to 504 >*/
  752.     if (xtf == 0.) {
  753.     goto L504;
  754.     }
  755. /*<       argtf=xtf >*/
  756.     argtf = xtf;
  757. /*<       if(ovtf.ne.0.0d0) argtf=argtf*dexp(vbc*ovtf) >*/
  758.     if (ovtf != 0.) {
  759.     argtf *= exp(vbc * ovtf);
  760.     }
  761. /*<       arg2=argtf >*/
  762.     arg2 = argtf;
  763. /*<       if(xjtf.eq.0.0d0) go to 503 >*/
  764.     if (xjtf == 0.) {
  765.     goto L503;
  766.     }
  767. /*<       temp=cbe/(cbe+xjtf) >*/
  768.     temp = cbe / (cbe + xjtf);
  769. /*<       argtf=argtf*temp*temp >*/
  770.     argtf = argtf * temp * temp;
  771. /*<       arg2=argtf*(3.0d0-temp-temp) >*/
  772.     arg2 = argtf * (3. - temp - temp);
  773. /*<   503 arg3=cbe*argtf*ovtf >*/
  774. L503:
  775.     arg3 = cbe * argtf * ovtf;
  776. /*<   504 cbe=cbe*(1.0d0+argtf)/qb >*/
  777. L504:
  778.     cbe = cbe * (argtf + 1.) / qb;
  779. /*<       gbe=(gbe*(1.0d0+arg2)-cbe*dqbdve)/qb >*/
  780.     gbe = (gbe * (arg2 + 1.) - cbe * dqbdve) / qb;
  781. /*<       geqcb=tf*(arg3-cbe*dqbdvc)/qb >*/
  782.     geqcb = tf * (arg3 - cbe * dqbdvc) / qb;
  783. /*<   505 if (vbe.ge.fcpe) go to 510 >*/
  784. L505:
  785.     if (vbe >= fcpe) {
  786.     goto L510;
  787.     }
  788. /*<       arg=1.0d0-vbe/pe >*/
  789.     arg = 1. - vbe / pe;
  790. /*<       sarg=dexp(-xme*dlog(arg)) >*/
  791.     sarg = exp(-xme * log(arg));
  792. /*<       qbe(lx0+loct)=tf*cbe+pe*czbe*(1.0d0-arg*sarg)/(1.0d0-xme) >*/
  793.     qbe[tabinf_1.lx0 + loct - 1] = tf * cbe + pe * czbe * (1. - arg * sarg) / 
  794.         (1. - xme);
  795. /*<       capbe=tf*gbe+czbe*sarg >*/
  796.     capbe = tf * gbe + czbe * sarg;
  797. /*<       go to 520 >*/
  798.     goto L520;
  799. /*<   510 f1=value(locm+47) >*/
  800. L510:
  801.     f1 = blank_1.value[locm + 46];
  802. /*<       f2=value(locm+48) >*/
  803.     f2 = blank_1.value[locm + 47];
  804. /*<       f3=value(locm+49) >*/
  805.     f3 = blank_1.value[locm + 48];
  806. /*<       czbef2=czbe/f2 >*/
  807.     czbef2 = czbe / f2;
  808. /*<       qbe(lx0+loct)=tf*cbe+czbe*f1+czbef2*(f3*(vbe-fcpe) >*/
  809. /*<      1   +(xme/(pe+pe))*(vbe*vbe-fcpe*fcpe)) >*/
  810.     qbe[tabinf_1.lx0 + loct - 1] = tf * cbe + czbe * f1 + czbef2 * (f3 * (vbe 
  811.         - fcpe) + xme / (pe + pe) * (vbe * vbe - fcpe * fcpe));
  812. /*<       capbe=tf*gbe+czbef2*(f3+xme*vbe/pe) >*/
  813.     capbe = tf * gbe + czbef2 * (f3 + xme * vbe / pe);
  814. /*<   520 fcpc=value(locm+50) >*/
  815. L520:
  816.     fcpc = blank_1.value[locm + 49];
  817. /*<       f1=value(locm+51) >*/
  818.     f1 = blank_1.value[locm + 50];
  819. /*<       f2=value(locm+52) >*/
  820.     f2 = blank_1.value[locm + 51];
  821. /*<       f3=value(locm+53) >*/
  822.     f3 = blank_1.value[locm + 52];
  823. /*<       if (vbc.ge.fcpc) go to 530 >*/
  824.     if (vbc >= fcpc) {
  825.     goto L530;
  826.     }
  827. /*<       arg=1.0d0-vbc/pc >*/
  828.     arg = 1. - vbc / pc;
  829. /*<       sarg=dexp(-xmc*dlog(arg)) >*/
  830.     sarg = exp(-xmc * log(arg));
  831. /*<       qbc(lx0+loct)=tr*cbc+pc*czbc*(1.0d0-arg*sarg)/(1.0d0-xmc) >*/
  832.     qbc[tabinf_1.lx0 + loct - 1] = tr * cbc + pc * czbc * (1. - arg * sarg) / 
  833.         (1. - xmc);
  834. /*<       capbc=tr*gbc+czbc*sarg >*/
  835.     capbc = tr * gbc + czbc * sarg;
  836. /*<       go to 540 >*/
  837.     goto L540;
  838. /*<   530 czbcf2=czbc/f2 >*/
  839. L530:
  840.     czbcf2 = czbc / f2;
  841. /*<       qbc(lx0+loct)=tr*cbc+czbc*f1+czbcf2*(f3*(vbc-fcpc) >*/
  842. /*<      1   +(xmc/(pc+pc))*(vbc*vbc-fcpc*fcpc)) >*/
  843.     qbc[tabinf_1.lx0 + loct - 1] = tr * cbc + czbc * f1 + czbcf2 * (f3 * (vbc 
  844.         - fcpc) + xmc / (pc + pc) * (vbc * vbc - fcpc * fcpc));
  845. /*<       capbc=tr*gbc+czbcf2*(f3+xmc*vbc/pc) >*/
  846.     capbc = tr * gbc + czbcf2 * (f3 + xmc * vbc / pc);
  847. /*<   540 if(vbx.ge.fcpc) go to 550 >*/
  848. L540:
  849.     if (vbx >= fcpc) {
  850.     goto L550;
  851.     }
  852. /*<       arg=1.0d0-vbx/pc >*/
  853.     arg = 1. - vbx / pc;
  854. /*<       sarg=dexp(-xmc*dlog(arg)) >*/
  855.     sarg = exp(-xmc * log(arg));
  856. /*<       qbx(lx0+loct)=pc*czbx*(1.0d0-arg*sarg)/(1.0d0-xmc) >*/
  857.     qbx[tabinf_1.lx0 + loct - 1] = pc * czbx * (1. - arg * sarg) / (1. - xmc);
  858.  
  859. /*<       capbx=czbx*sarg >*/
  860.     capbx = czbx * sarg;
  861. /*<       go to 560 >*/
  862.     goto L560;
  863. /*<   550 czbxf2=czbx/f2 >*/
  864. L550:
  865.     czbxf2 = czbx / f2;
  866. /*<       qbx(lx0+loct)=czbx*f1+czbxf2*(f3*(vbx-fcpc)+(xmc/(pc+pc))* >*/
  867. /*<      1   (vbx*vbx-fcpc*fcpc)) >*/
  868.     qbx[tabinf_1.lx0 + loct - 1] = czbx * f1 + czbxf2 * (f3 * (vbx - fcpc) + 
  869.         xmc / (pc + pc) * (vbx * vbx - fcpc * fcpc));
  870. /*<       capbx=czbxf2*(f3+xmc*vbx/pc) >*/
  871.     capbx = czbxf2 * (f3 + xmc * vbx / pc);
  872. /*<   560 if(vcs.ge.0.0d0) go to 570 >*/
  873. L560:
  874.     if (vcs >= 0.) {
  875.     goto L570;
  876.     }
  877. /*<       arg=1.0d0-vcs/ps >*/
  878.     arg = 1. - vcs / ps;
  879. /*<       sarg=dexp(-xms*dlog(arg)) >*/
  880.     sarg = exp(-xms * log(arg));
  881. /*<       qcs(lx0+loct)=ps*czcs*(1.0d0-arg*sarg)/(1.0d0-xms) >*/
  882.     qcs[tabinf_1.lx0 + loct - 1] = ps * czcs * (1. - arg * sarg) / (1. - xms);
  883.  
  884. /*<       capcs=czcs*sarg >*/
  885.     capcs = czcs * sarg;
  886. /*<       go to 580 >*/
  887.     goto L580;
  888. /*<   570 qcs(lx0+loct)=vcs*czcs*(1.0d0+xms*vcs/(2.0d0*ps)) >*/
  889. L570:
  890.     qcs[tabinf_1.lx0 + loct - 1] = vcs * czcs * (xms * vcs / (ps * 2.) + 1.);
  891. /*<       capcs=czcs*(1.0d0+xms*vcs/ps) >*/
  892.     capcs = czcs * (xms * vcs / ps + 1.);
  893.  
  894. /*  store small-signal parameters */
  895.  
  896. /*<   580 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700 >*/
  897. L580:
  898.     if (status_1.mode == 1 && status_1.modedc == 2 && status_1.nosolv != 0) {
  899.     goto L700;
  900.     }
  901. /*<       if (initf.ne.4) go to 600 >*/
  902.     if (status_1.initf != 4) {
  903.     goto L600;
  904.     }
  905. /*<       value(lx0+loct+9)=capbe >*/
  906.     blank_1.value[tabinf_1.lx0 + loct + 8] = capbe;
  907. /*<       value(lx0+loct+11)=capbc >*/
  908.     blank_1.value[tabinf_1.lx0 + loct + 10] = capbc;
  909. /*<       value(lx0+loct+13)=capcs >*/
  910.     blank_1.value[tabinf_1.lx0 + loct + 12] = capcs;
  911. /*<       value(lx0+loct+15)=capbx >*/
  912.     blank_1.value[tabinf_1.lx0 + loct + 14] = capbx;
  913. /*<       value(lx0+loct+17)=geqcb >*/
  914.     blank_1.value[tabinf_1.lx0 + loct + 16] = geqcb;
  915. /*<       go to 1000 >*/
  916.     goto L1000;
  917.  
  918. /*  transient analysis */
  919.  
  920. /*<   600 if (initf.ne.5) go to 610 >*/
  921. L600:
  922.     if (status_1.initf != 5) {
  923.     goto L610;
  924.     }
  925. /*<       qbe(lx1+loct)=qbe(lx0+loct) >*/
  926.     qbe[tabinf_1.lx1 + loct - 1] = qbe[tabinf_1.lx0 + loct - 1];
  927. /*<       qbc(lx1+loct)=qbc(lx0+loct) >*/
  928.     qbc[tabinf_1.lx1 + loct - 1] = qbc[tabinf_1.lx0 + loct - 1];
  929. /*<       qbx(lx1+loct)=qbx(lx0+loct) >*/
  930.     qbx[tabinf_1.lx1 + loct - 1] = qbx[tabinf_1.lx0 + loct - 1];
  931. /*<       qcs(lx1+loct)=qcs(lx0+loct) >*/
  932.     qcs[tabinf_1.lx1 + loct - 1] = qcs[tabinf_1.lx0 + loct - 1];
  933. /*<   610 call intgr8(geq,ceq,capbe,loct+8) >*/
  934. L610:
  935.     i_1 = loct + 8;
  936.     intgr8_(&geq, &ceq, &capbe, &i_1);
  937. /*<       geqcb=geqcb*ag(1) >*/
  938.     geqcb *= status_1.ag[0];
  939. /*<       gpi=gpi+geq >*/
  940.     gpi += geq;
  941. /*<       cb=cb+cqbe(lx0+loct) >*/
  942.     cb += cqbe[tabinf_1.lx0 + loct - 1];
  943. /*<       call intgr8(geq,ceq,capbc,loct+10) >*/
  944.     i_1 = loct + 10;
  945.     intgr8_(&geq, &ceq, &capbc, &i_1);
  946. /*<       gmu=gmu+geq >*/
  947.     gmu += geq;
  948. /*<       cb=cb+cqbc(lx0+loct) >*/
  949.     cb += cqbc[tabinf_1.lx0 + loct - 1];
  950. /*<       cc=cc-cqbc(lx0+loct) >*/
  951.     cc -= cqbc[tabinf_1.lx0 + loct - 1];
  952. /*<       if (initf.ne.5) go to 700 >*/
  953.     if (status_1.initf != 5) {
  954.     goto L700;
  955.     }
  956. /*<       cqbe(lx1+loct)=cqbe(lx0+loct) >*/
  957.     cqbe[tabinf_1.lx1 + loct - 1] = cqbe[tabinf_1.lx0 + loct - 1];
  958. /*<       cqbc(lx1+loct)=cqbc(lx0+loct) >*/
  959.     cqbc[tabinf_1.lx1 + loct - 1] = cqbc[tabinf_1.lx0 + loct - 1];
  960.  
  961. /*  check convergence */
  962.  
  963. /*<   700 if (initf.ne.3) go to 710 >*/
  964. L700:
  965.     if (status_1.initf != 3) {
  966.     goto L710;
  967.     }
  968. /*<       if (ioff.eq.0) go to 710 >*/
  969.     if (ioff == 0) {
  970.     goto L710;
  971.     }
  972. /*<       go to 750 >*/
  973.     goto L750;
  974. /*<   710 if (icheck.eq.1) go to 720 >*/
  975. L710:
  976.     if (icheck == 1) {
  977.     goto L720;
  978.     }
  979. /*<       tol=reltol*dmax1(dabs(cchat),dabs(cc))+abstol >*/
  980. /* Computing MAX */
  981.     d_1 = abs(cchat), d_2 = abs(cc);
  982.     tol = knstnt_1.reltol * max(d_2,d_1) + knstnt_1.abstol;
  983. /*<       if (dabs(cchat-cc).gt.tol) go to 720 >*/
  984.     if ((d_1 = cchat - cc, abs(d_1)) > tol) {
  985.     goto L720;
  986.     }
  987. /*<       tol=reltol*dmax1(dabs(cbhat),dabs(cb))+abstol >*/
  988. /* Computing MAX */
  989.     d_1 = abs(cbhat), d_2 = abs(cb);
  990.     tol = knstnt_1.reltol * max(d_2,d_1) + knstnt_1.abstol;
  991. /*<       if (dabs(cbhat-cb).le.tol) go to 750 >*/
  992.     if ((d_1 = cbhat - cb, abs(d_1)) <= tol) {
  993.     goto L750;
  994.     }
  995. /*<   720 noncon=noncon+1 >*/
  996. L720:
  997.     ++status_1.noncon;
  998. /*<   750 vbeo(lx0+loct)=vbe >*/
  999. L750:
  1000.     vbeo[tabinf_1.lx0 + loct - 1] = vbe;
  1001. /*<       vbco(lx0+loct)=vbc >*/
  1002.     vbco[tabinf_1.lx0 + loct - 1] = vbc;
  1003. /*<       cco(lx0+loct)=cc >*/
  1004.     cco[tabinf_1.lx0 + loct - 1] = cc;
  1005. /*<       cbo(lx0+loct)=cb >*/
  1006.     cbo[tabinf_1.lx0 + loct - 1] = cb;
  1007. /*<       gpio(lx0+loct)=gpi >*/
  1008.     gpio[tabinf_1.lx0 + loct - 1] = gpi;
  1009. /*<       gmuo(lx0+loct)=gmu >*/
  1010.     gmuo[tabinf_1.lx0 + loct - 1] = gmu;
  1011. /*<       gmo(lx0+loct)=gm >*/
  1012.     gmo[tabinf_1.lx0 + loct - 1] = gm;
  1013. /*<       goo(lx0+loct)=go >*/
  1014.     goo[tabinf_1.lx0 + loct - 1] = go;
  1015. /*<       gxo(lx0+loct)=gx >*/
  1016.     gxo[tabinf_1.lx0 + loct - 1] = gx;
  1017. /*<       geqcbo(lx0+loct)=geqcb >*/
  1018.     geqcbo[tabinf_1.lx0 + loct - 1] = geqcb;
  1019. /*<       if (mode.eq.1) go to 900 >*/
  1020.     if (status_1.mode == 1) {
  1021.     goto L900;
  1022.     }
  1023.  
  1024. /*     charge storage for c-s and b-x junctions */
  1025.  
  1026. /*<   800 call intgr8(gccs,ceq,capcs,loct+12) >*/
  1027. L800:
  1028.     i_1 = loct + 12;
  1029.     intgr8_(&gccs, &ceq, &capcs, &i_1);
  1030. /*<       ceqcs=type*(cqcs(lx0+loct)-vcs*gccs) >*/
  1031.     ceqcs = type * (cqcs[tabinf_1.lx0 + loct - 1] - vcs * gccs);
  1032. /*<       call intgr8(geqbx,ceq,capbx,loct+14) >*/
  1033.     i_1 = loct + 14;
  1034.     intgr8_(&geqbx, &ceq, &capbx, &i_1);
  1035. /*<       ceqbx=type*(cqbx(lx0+loct)-vbx*geqbx) >*/
  1036.     ceqbx = type * (cqbx[tabinf_1.lx0 + loct - 1] - vbx * geqbx);
  1037. /*<       if (initf.ne.5) go to 900 >*/
  1038.     if (status_1.initf != 5) {
  1039.     goto L900;
  1040.     }
  1041. /*<       cqbx(lx1+loct)=cqbx(lx0+loct) >*/
  1042.     cqbx[tabinf_1.lx1 + loct - 1] = cqbx[tabinf_1.lx0 + loct - 1];
  1043. /*<       cqcs(lx1+loct)=cqcs(lx0+loct) >*/
  1044.     cqcs[tabinf_1.lx1 + loct - 1] = cqcs[tabinf_1.lx0 + loct - 1];
  1045.  
  1046. /*  load current excitation vector */
  1047.  
  1048. /*<   900 ceqbe=type*(cc+cb-vbe*(gm+go+gpi)+vbc*(go-geqcb)) >*/
  1049. L900:
  1050.     ceqbe = type * (cc + cb - vbe * (gm + go + gpi) + vbc * (go - geqcb));
  1051. /*<       ceqbc=type*(-cc+vbe*(gm+go)-vbc*(gmu+go)) >*/
  1052.     ceqbc = type * (-cc + vbe * (gm + go) - vbc * (gmu + go));
  1053. /*<       value(lvn+node2)=value(lvn+node2)-ceqbx >*/
  1054.     blank_1.value[tabinf_1.lvn + node2 - 1] -= ceqbx;
  1055. /*<       value(lvn+node4)=value(lvn+node4)+ceqcs+ceqbx+ceqbc >*/
  1056.     blank_1.value[tabinf_1.lvn + node4 - 1] = blank_1.value[tabinf_1.lvn + 
  1057.         node4 - 1] + ceqcs + ceqbx + ceqbc;
  1058. /*<       value(lvn+node5)=value(lvn+node5)-ceqbe-ceqbc >*/
  1059.     blank_1.value[tabinf_1.lvn + node5 - 1] = blank_1.value[tabinf_1.lvn + 
  1060.         node5 - 1] - ceqbe - ceqbc;
  1061. /*<       value(lvn+node6)=value(lvn+node6)+ceqbe >*/
  1062.     blank_1.value[tabinf_1.lvn + node6 - 1] += ceqbe;
  1063. /*<       value(lvn+node7)=value(lvn+node7)-ceqcs >*/
  1064.     blank_1.value[tabinf_1.lvn + node7 - 1] -= ceqcs;
  1065.  
  1066. /*  load y matrix */
  1067.  
  1068. /*<       locy=lvn+nodplc(loc+24) >*/
  1069.     locy = tabinf_1.lvn + nodplc[loc + 23];
  1070. /*<       value(locy)=value(locy)+gcpr >*/
  1071.     blank_1.value[locy - 1] += gcpr;
  1072. /*<       locy=lvn+nodplc(loc+25) >*/
  1073.     locy = tabinf_1.lvn + nodplc[loc + 24];
  1074. /*<       value(locy)=value(locy)+gx+geqbx >*/
  1075.     blank_1.value[locy - 1] = blank_1.value[locy - 1] + gx + geqbx;
  1076. /*<       locy=lvn+nodplc(loc+26) >*/
  1077.     locy = tabinf_1.lvn + nodplc[loc + 25];
  1078. /*<       value(locy)=value(locy)+gepr >*/
  1079.     blank_1.value[locy - 1] += gepr;
  1080. /*<       locy=lvn+nodplc(loc+27) >*/
  1081.     locy = tabinf_1.lvn + nodplc[loc + 26];
  1082. /*<       value(locy)=value(locy)+gmu+go+gcpr+gccs+geqbx >*/
  1083.     blank_1.value[locy - 1] = blank_1.value[locy - 1] + gmu + go + gcpr + 
  1084.         gccs + geqbx;
  1085. /*<       locy=lvn+nodplc(loc+28) >*/
  1086.     locy = tabinf_1.lvn + nodplc[loc + 27];
  1087. /*<       value(locy)=value(locy)+gx  +gpi+gmu+geqcb >*/
  1088.     blank_1.value[locy - 1] = blank_1.value[locy - 1] + gx + gpi + gmu + 
  1089.         geqcb;
  1090. /*<       locy=lvn+nodplc(loc+29) >*/
  1091.     locy = tabinf_1.lvn + nodplc[loc + 28];
  1092. /*<       value(locy)=value(locy)+gpi+gepr+gm+go >*/
  1093.     blank_1.value[locy - 1] = blank_1.value[locy - 1] + gpi + gepr + gm + go;
  1094. /*<       locy=lvn+nodplc(loc+10) >*/
  1095.     locy = tabinf_1.lvn + nodplc[loc + 9];
  1096. /*<       value(locy)=value(locy)-gcpr >*/
  1097.     blank_1.value[locy - 1] -= gcpr;
  1098. /*<       locy=lvn+nodplc(loc+11) >*/
  1099.     locy = tabinf_1.lvn + nodplc[loc + 10];
  1100. /*<       value(locy)=value(locy)-gx >*/
  1101.     blank_1.value[locy - 1] -= gx;
  1102. /*<       locy=lvn+nodplc(loc+12) >*/
  1103.     locy = tabinf_1.lvn + nodplc[loc + 11];
  1104. /*<       value(locy)=value(locy)-gepr >*/
  1105.     blank_1.value[locy - 1] -= gepr;
  1106. /*<       locy=lvn+nodplc(loc+13) >*/
  1107.     locy = tabinf_1.lvn + nodplc[loc + 12];
  1108. /*<       value(locy)=value(locy)-gcpr >*/
  1109.     blank_1.value[locy - 1] -= gcpr;
  1110. /*<       locy=lvn+nodplc(loc+14) >*/
  1111.     locy = tabinf_1.lvn + nodplc[loc + 13];
  1112. /*<       value(locy)=value(locy)-gmu+gm >*/
  1113.     blank_1.value[locy - 1] = blank_1.value[locy - 1] - gmu + gm;
  1114. /*<       locy=lvn+nodplc(loc+15) >*/
  1115.     locy = tabinf_1.lvn + nodplc[loc + 14];
  1116. /*<       value(locy)=value(locy)-gm-go >*/
  1117.     blank_1.value[locy - 1] = blank_1.value[locy - 1] - gm - go;
  1118. /*<       locy=lvn+nodplc(loc+16) >*/
  1119.     locy = tabinf_1.lvn + nodplc[loc + 15];
  1120. /*<       value(locy)=value(locy)-gx >*/
  1121.     blank_1.value[locy - 1] -= gx;
  1122. /*<       locy=lvn+nodplc(loc+17) >*/
  1123.     locy = tabinf_1.lvn + nodplc[loc + 16];
  1124. /*<       value(locy)=value(locy)-gmu-geqcb >*/
  1125.     blank_1.value[locy - 1] = blank_1.value[locy - 1] - gmu - geqcb;
  1126. /*<       locy=lvn+nodplc(loc+18) >*/
  1127.     locy = tabinf_1.lvn + nodplc[loc + 17];
  1128. /*<       value(locy)=value(locy)-gpi >*/
  1129.     blank_1.value[locy - 1] -= gpi;
  1130. /*<       locy=lvn+nodplc(loc+19) >*/
  1131.     locy = tabinf_1.lvn + nodplc[loc + 18];
  1132. /*<       value(locy)=value(locy)-gepr >*/
  1133.     blank_1.value[locy - 1] -= gepr;
  1134. /*<       locy=lvn+nodplc(loc+20) >*/
  1135.     locy = tabinf_1.lvn + nodplc[loc + 19];
  1136. /*<       value(locy)=value(locy)-go+geqcb >*/
  1137.     blank_1.value[locy - 1] = blank_1.value[locy - 1] - go + geqcb;
  1138. /*<       locy=lvn+nodplc(loc+21) >*/
  1139.     locy = tabinf_1.lvn + nodplc[loc + 20];
  1140. /*<       value(locy)=value(locy)-gpi-gm-geqcb >*/
  1141.     blank_1.value[locy - 1] = blank_1.value[locy - 1] - gpi - gm - geqcb;
  1142. /*<       locy=lvn+nodplc(loc+31) >*/
  1143.     locy = tabinf_1.lvn + nodplc[loc + 30];
  1144. /*<       value(locy)=value(locy)+gccs >*/
  1145.     blank_1.value[locy - 1] += gccs;
  1146. /*<       locy=lvn+nodplc(loc+32) >*/
  1147.     locy = tabinf_1.lvn + nodplc[loc + 31];
  1148. /*<       value(locy)=value(locy)-gccs >*/
  1149.     blank_1.value[locy - 1] -= gccs;
  1150. /*<       locy=lvn+nodplc(loc+33) >*/
  1151.     locy = tabinf_1.lvn + nodplc[loc + 32];
  1152. /*<       value(locy)=value(locy)-gccs >*/
  1153.     blank_1.value[locy - 1] -= gccs;
  1154. /*<       locy=lvn+nodplc(loc+34) >*/
  1155.     locy = tabinf_1.lvn + nodplc[loc + 33];
  1156. /*<       value(locy)=value(locy)-geqbx >*/
  1157.     blank_1.value[locy - 1] -= geqbx;
  1158. /*<       locy=lvn+nodplc(loc+35) >*/
  1159.     locy = tabinf_1.lvn + nodplc[loc + 34];
  1160. /*<       value(locy)=value(locy)-geqbx >*/
  1161.     blank_1.value[locy - 1] -= geqbx;
  1162. /*<  1000 loc=nodplc(loc) >*/
  1163. L1000:
  1164.     loc = nodplc[loc - 1];
  1165. /*<       go to 10 >*/
  1166.     goto L10;
  1167. /*<       end >*/
  1168. } /* bjt_ */
  1169.  
  1170. #undef gxo
  1171. #undef qbx
  1172. #undef qcs
  1173. #undef goo
  1174. #undef gmo
  1175. #undef qbe
  1176. #undef qbc
  1177. #undef cco
  1178. #undef cbo
  1179. #undef cvalue
  1180. #undef nodplc
  1181. #undef geqcbo
  1182. #undef cexbc
  1183. #undef gmuo
  1184. #undef cqbx
  1185. #undef gpio
  1186. #undef cqcs
  1187. #undef vbeo
  1188. #undef vbco
  1189. #undef cqbe
  1190. #undef cqbc
  1191.  
  1192.  
  1193.